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We  present  a  numerical  investigation  of  the  second-order  nonlinear  optical  properties  of  metal-based 
metamaterial  nanoresonators.  The  nonlinear  optical  response  of  the  metal  is  described  by  a  hydrodynamic 
model,  with  the  effects  of  electron  pressure  in  the  electron  gas  also  taken  into  account.  We  show  that  as  the 
pressure  term  tends  to  zero  the  amount  of  converted  second-harmonic  field  tends  to  an  asymptotic  value.  In  this 
limit  it  becomes  possible  to  rewrite  the  nonlinear  surface  contributions  as  functions  of  the  value  of  the  polarization 
vector  inside  the  bulk  region.  Nonlocality  thus  can  be  incorporated  into  numerical  simulations  without  actually 
utilizing  the  nonlocal  equation  of  motion  or  solving  for  the  rapidly  varying  fields  that  occur  near  the  metal 
surface.  We  use  our  model  to  investigate  the  second-harmonic  generation  process  with  three-dimensional  gold 
nanoparticle  arrays  and  show  that  nanocrescents  can  easily  attain  conversion  efficiencies  of  -6.0  x  10-8  for 
pumping  peak  intensities  of  a  few  tens  of  MW/cm2. 
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I.  INTRODUCTION 

Metamaterials  are  artificially  structured  media  whose  col¬ 
lective  electromagnetic  properties  can  differ  markedly  from 
those  displayed  by  the  individual  components.  The  advent 
of  metamaterials  has  introduced  a  paradigm  shift  in  the 
concept  of  what  constitutes  a  material  and  what  determines  its 
fundamental  electromagnetic  and  other  physical  properties. 
The  possibility  of  tailoring  magnetic  and  electric  material 
responses  at  will,  rather  than  merely  exploiting  materials 
provided  by  nature,  has  enabled  the  pursuit  of  unusual 
electromagnetic  properties,  such  as  invisibility1  and  negative 
refractive  index.2 

Without  a  doubt,  metamaterials  constitute  novel  and  in¬ 
teresting  linear  media  that  have  provided  a  venue  to  explore 
otherwise  inaccessible  concepts.  Yet,  artificial  materials  in  the 
context  of  nonlinear  optical  media  may  offer  even  greater 
opportunities,  since  the  inherent  local  field  inhomogeneities 
associated  with  metamaterial  inclusions  can  translate  into  large 
enhancement  of  the  local  field  and  significant  lowering  of  the 
nonlinear  thresholds. 

Nonlinear  metamaterial  composites  have  been  analyzed 
both  theoretically  and  experimentally  in  the  microwave  and 
near-infrared  regions  of  the  spectrum.3-5  For  such  materials, 
one  may  apply  an  analytical  homogenization  model  that 
leads  to  closed-form  expressions  for  the  effective  nonlinear 
susceptibilities,6  placing  the  design  of  nonlinear  metamaterials 
on  an  equal  footing  with  linear  metamaterials.  In  addition, 
numerical  retrieval  approaches  have  been  developed  that 
allow  precise  values  to  be  ascribed  to  the  effective  nonlinear 
susceptibilities  for  a  composite  metamaterial  comprising 
both  structured  inclusions  as  well  as  embedded  nonlinear 
elements.7,8 

The  most  common  metamaterial  designs  have  made  use  of 
conducting  metal  inclusions  that  function  as  subwavelength 
electrical  circuits.  At  low  frequencies,  metals  behave  as 
conductors,  and,  thus,  metal  inclusions  can  be  conceptually 
divided  into  inductive,  capacitive,  and  resistive  regions,  with  a 
resonance  frequency  determined  by  the  inductance  and  capac¬ 


itance  in  the  usual  manner.  If  any  of  these  circuit  parameters 
are  made  nonlinear — through  the  introduction  of  an  actual 
lumped  component  such  as  a  varactor  diode  or  an  inherently 
nonlinear  crystal  embedded  into  the  capacitive  region  of  the 
inclusion — the  effective  circuit  then  produces  a  nonlinear 
response  to  the  driving  electric  or  magnetic  field  (depending  on 
how  the  wave  couples  to  the  inclusion).  Applying  analytical 
or  retrieval  techniques  on  such  a  structure  yields  adequate 
effective  medium  values  for  the  linear  constitutive  tensor 
elements,  as  well  as  for  the  tensor  elements  of  the  nonlinear 
susceptibility  terms.  A  nonlinear  metamaterial  is  advantageous 
in  that  nonlinear  thresholds  may  be  significantly  lowered,  often 
irrespective  of  the  linear  properties  of  the  composite.  That  is, 
the  linear  and  nonlinear  responses  of  a  nonlinear  metamaterial 
can  be  designed  with  considerable  independence,  allowing 
customized  anisotropy  and  other  properties  that  may  be  used 
to  improve  the  efficiency  of  nonlinear  applications. 

Nonlinear  metamaterials  based  on  metals  are,  thus,  ap¬ 
pealing  for  potential  use  in  nonlinear  optical  applications 
at  infrared,  visible,  and  ultraviolet  wavelengths.  However,  a 
simple  scaling  of  the  geometrical  parameters  of  metamaterial 
inclusions  from  microwave  to  visible  frequencies,  for  example, 
is  usually  not  sufficient  to  create  a  viable  nonlinear  optical 
metamaterial.  At  frequencies  above  a  few  terahertz,  metal 
response  changes  from  conductor-like  to  dielectric-like,  with 
considerable  absorption  occurring  as  the  electromagnetic 
radiation  extends  further  and  further  into  the  metal.  Though 
the  electrical  circuit  properties  of  metal  inclusions  persist  at 
visible  wavelengths,  the  inertia  of  the  charge  carriers  becomes 
an  increasingly  dominant  contribution  to  the  inductance, 
and  absorption  increases  significantly.9  The  design  of  metal 
optical  metamaterials — where  optical  is  roughly  defined  as 
being  above  10-12  THz — requires  a  different  approach,  with 
different  applications  being  favored. 

In  addition  to  field  enhancement  being  a  useful  mecha¬ 
nism  in  optical  metamaterials,  the  intrinsic  nonlinearity  of 
metals  makes  metal-based  nonlinear  optical  metamaterials  an 
interesting  possibility.  For  example,  values  of  / (3)  for  metals 
at  optical  wavelengths  are  competitive  with  most  materials. 
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X(3)  of  gold  is  5.4  x  10“ 19  m2/V2,  somewhat  lower  than 
semiconductors  such  as  silicon  (x(3)  =  2.0  x  10-18  m2/V2) 
but  larger  than  other  common  crystals  such  as  diamond 
(X(3)  =  1.8  x  10“21  m2/V2)  or  A1203  (x(3)  =  2.2  x  10“22 
m2/V2).10  Moreover,  though  metals  are  centrosymmetric  and 
do  not  possess  an  inherent  x(2)  nonlinearity,  the  surface  of  a 
metal  can  break  the  symmetry  and  provide  a  mechanism  for 
an  effective  x (2)  nonlinearity.  The  homogenized  x (2)  nonlinear 
response  of  a  metal,  thus,  becomes  highly  dependent  on  its 
geometry,  making  it  inherently  a  metamaterial  construct.  In 
all  cases,  the  large  absorption  associated  with  metals  suggests 
nanostructuring  as  a  means  of  optimizing  the  nonlinear  effects 
while  minimizing  propagation  within  the  metal.  As  we  show 
here,  nonlinear  metal-based  metamaterials  can  accomplish 
these  goals. 

The  origin  of  the  nonlinearity  in  metals  arises  from  the 
response  of  both  bound  and  free  electrons.  In  particular, 
for  the  visible/near-IR  part  of  the  spectrum,  the  nonlinear 
response  of  thick  metal  layers  may  be  attributed  mostly  to  free 
electrons,  with  Lorentz  (magnetic)  and  quadrupolar  contribu¬ 
tions  from  bound  charges  becoming  increasingly  important 
at  shorter  wavelengths.  The  interest  in  nonlinear  interactions 
at  metal  surfaces  dates  back  to  the  beginning  of  nonlinear 
optics.11-17  Optical  second-harmonic  generation  from  a  silver 
slab  was  first  observed  in  1965. 18,19  Since  then,  numerous 
experimental  studies  have  been  published,20-24  with  a  variety 
of  phenomenological25-33  and  microscopic  approaches34-37 
developed  to  analyze  the  second-harmonic  response  of  metal 
structures. 

Recent  research  in  plasmonic  phenomena  has  stimulated 
renewed  interest  in  the  nonlinear  optical  properties  of  metals. 
In  particular,  a  steady  stream  of  works  concerning  harmonic 
generation  from  metallic  nanostructures,  such  as  hole  ar¬ 
rays  in  a  metallic  substrate,38-43  metallodielectric  multilayer 
structures,44,45  Au  nanoantennas, 46,47  and  periodic  nanostruc- 
tured  metal  films,48,49  have  been  published  recently.  Second- 
harmonic  generation  has  also  been  observed  experimentally 
from  ordered  arrays  of  split-ring  resonators,50,51  as  well  as 
from  a  variety  of  single  nanoparticles52-56  and  nanoparticle 
arrays  with  varying  geometries.57-69  On  the  theoretical  side, 
Zeng  et  al.  readapted  the  free-electron  theory  of  second- 
harmonic  generation  to  arbitrary  shaped  metal  nanoparticles.70 
Scalora  et  al.  have  considered  a  dynamic  description  based 
on  the  hydrodynamic  model  and  taken  into  account  bound 
electron  contributions  to  the  second-  and  third-harmonic 
generation  in  the  ultrashort  pulse  regime.71 

It  should  be  emphasized  that  metals  contain  both  volume 
and  surface  nonlinear  contributions,28  and  it  is  generally 
not  possible  to  entirely  separate  the  two.  The  separation 
becomes  even  more  ambiguous  for  subwavelength  nanolayers 
or  nanoparticles,  where  the  field  can  penetrate  the  metal 
itself  so  the  nonlinear  response  may  be  distributed  inside 
the  volume.  Under  these  circumstances  the  shape  of  the 
nanoparticle  acquires  a  more  crucial  role  in  determining 
the  nonlinear  response  of  the  metal.  Moreover,  nonlinear 
surface  contributions  are  strictly  related  to  the  response 
of  the  electrons  within  the  Thomas -Fermi  screening  length 
(A.tf  ~  1  A  for  gold)  from  the  surface,  which  introduces  a 
crucial  microscopic  scale  to  the  problem.  Thus,  solving  for 
the  optical  fields  in  the  macroscopic  realm,  which  vary  on 


the  order  of  several  hundred  nanometers,  can  rapidly  result  in 
considerable  numerical  complexity. 

In  this  work  we  perform  an  analysis  of  second-harmonic 
generation  in  plasmonic  systems  of  arbitrary  shape.  We, 
first,  summarize  the  hydrodynamic  model  that  describes  the 
nonlinear  optical  response  of  the  metal,  including  the  effects 
of  quantum  pressure  associated  with  the  electron  gas.  A  numer¬ 
ical  analysis  of  the  resulting  nonlocal  and  nonlinear  problem 
is  performed,  with  particular  emphasis  on  surface  effects.  In 
particular,  we  numerically  investigate  the  nonlinear  response 
at  the  metal  surface  as  the  electron  pressure  tends  to  zero, 
thus  establishing  an  unambiguous  bridge  to  the  free-electron 
limit ,  in  which  the  pressure  is  completely  neglected.  In  a 
recent  work72  we  provided  a  formula  expressing  the  nonlinear 
surface  contribution  in  terms  of  the  polarization  vector  in  the 
bulk  regions.  Here  we  show  a  detailed  comparison  with  the 
full  resolution  of  nonlocal  and  nonlinear  equations.  We  then 
numerically  investigate  second-harmonic  generation  arising 
from  a  variety  of  three-dimensional  nanoparticle  arrays  and 
show  the  impact  of  specific  geometries,  such  as  nanocrescents, 
on  the  enhancement  of  conversion  efficiencies. 


II.  MODEL 


A  description  for  the  polarization  inside  the  metal  is  given 
by  the  hydrodynamic  model.  Here  we  recall  the  most  salient 
elements  of  the  theory  and  refer  the  reader  to  Refs.  73,  35, 
and  71  for  further  details. 

The  electron  fluid  density,  n(r,t ),  and  the  current  density, 
J  =  env,  satisfy  Euler’s  equation, 


m^n 


3v 

Yt 


+  (v  •  V)v+  yy 


—  enE  +  en\  x  B  —  Wp, 


(1) 


along  with  the  continuity  equation, 


V  •  J  =  —eh, 


(2) 


where  the  dot  represents  the  partial  derivative  with  respect 
to  time,  m*  is  the  effective  electron  mass,  y  is  the  electron 
collision  rate,  v  is  the  electron  velocity  field,  and  p  is  the 
electron  pressure  which,  for  a  three-dimensional  gas,  is  given 

by71,74 


p(r,t)  =  po 


n(r,t) 


5/3 


n  o 


(3) 


where  po  ~  (EF  is  the  Fermi  energy  and  no  is  the 
equilibrium  charge  density).  Combining  Eqs.  (2),  (1),  and  (3), 
and  taking  into  account  P  =  J,  one  finds  that  the  free-electron 
polarization  P  satisfies  the  following  equation71 : 


p  +  yp=  — E-—  E(V-P)+—  PxB 
m*  m*  m* 

[(V  •  P)P  +  (P  •  V)P]  +  (V  •  P) 

noe  3  m* 

10  Ef 

-  - - “7  (V  •  P)  V  (V  •  P) .  (4) 

9  enoml 

In  calculating  Eq.  (4),  only  first-  and  second-order  terms  have 
been  retained.  Equation  (4)  summarizes  the  electron  fluid 
response  under  the  influence  of  the  electromagnetic  field.  In 
addition  to  the  linear  Drude  response,  we  have  the  nonlinear 
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Coulomb  term  (referred  to  as  quadrupole-like  term  by  virtue 
of  its  form),  proportional  to  E  (V  •  P),  the  magnetic  Lorentz 
force  contribution,  P  x  B,  the  convective  terms  (V  •  P)P  and 
(P  •  V)P,  the  linear  nonlocal  pressure  term  proportional  to 
V(V  •  P),  and  the  nonlinear  pressure  term  (V  •  P)V(V  •  P). 

In  order  to  calculate  the  expressions  for  the  fundamental 
and  second-harmonic  polarization,  we  expand  all  fields  in  a 
perturbative  manner, 

E(r,t)  =  Ei(r)e-i0>‘  +  E2(r)e-2io)'  +  ■■■, 

H(r,0  =  H1(r)e-,w?  +  n2(r)e~2io>t  +  •  •  • ,  (5) 

P(r,0  =  Pi(r)e-ia"  +  P2(r)e~2ia,t  +  ■■■, 

where  the  time  dependence  has  been  explicitly  indicated. 
Considering  terms  up  to  the  second  order,  we  obtain  the 
following  set  of  equations: 

2 

£2V(V  •  Pi)  +  («2  +  «o>)/)Pi  =  — —  Ei,  (6a) 

m*e 

2 

P2  V(V  •  P2)  +  (4 «2  +  2iojy )P2  =  -  — E2  +  SNL,  (6b) 

where  ft2  =  (5/3)EF/ml  and  the  nonlinear  source,  Snl,  is 
given  by 

e  icoe 

Snl  =  —Ei  (V  •  Pi)  H - -Pi  x  Bi 

m*  m* 

[(V-POPj  +  o*!  •  V)Pj] 
n0e 

10  Ef 

+  --^7(V.Pi)V(V.P1).  (7) 

9  enom* 

The  quantity  ft  is  proportional  to  the  Fermi  velocity  vF 
and,  consequently,  to  the  Thomas-Fermi  screening  length 
A.tf  =  vp/eop,  with  cov  being  the  plasma  frequency.  Equations 
(6)  with  Eq.  (7)  describe  the  fundamental  and  second- 
harmonic  polarization  vectors,  respectively,  and  hold  under 
the  assumption  that  the  fundamental  fields  are  not  affected  by 
the  generated  harmonic  (nondepleted  pump  approximation). 
However,  should  one  desire  to  calculate  the  self-consistent 
field  interaction,  the  necessary  terms  may  be  taken  into  account 
by  modifying  Eq.  (6a). 

Equations  (6)  are  solved  along  with  Maxwell’s  equations 
that,  in  the  framework  of  harmonic  propagation,  read75 

V  x  V  x  Ei  —  k2E\  =  iiqco2 Pi,  (8a) 

V  x  V  x  E2  —  k^E 2  =  iM)4co2F2,  (8b) 

where  k\  =  co/c  and  k2  =  2co/c,  with  c  the  light  velocity  in 
vacuum. 

The  polarization  defined  by  Eqs.  (6)  contain  linear,  nonlocal 
contributions  arising  from  the  electron  pressure  of  the  form 
/32V(V  •  P).  This  term  has  been  shown  to  introduce  significant 
modifications  in  the  optical  properties  of  subnanometer-gap 
metal-based  systems.76  It  has  also  been  predicted  to  be 
responsible  for  an  unusual,  resonant-like  phenomenon  for 
nano  wires  only  a  few  nanometers  in  diameter. 77-79  More 
generally,  this  term  takes  into  account  variations  of  the  electric 
field  occurring  in  a  region  of  the  order  of  A,TF  in  the  vicinity  of 
the  metal  surface,  where  electron-electron  interactions  become 
significant.  The  presence  of  spatial  derivatives  (nonlocality) 
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in  the  description  of  the  polarization  vector  requires  the 
specification  of  additional  boundary  conditions  to  solve  the 
electromagnetic  boundary  value  problem.80,81 

The  choice  of  additional  boundary  conditions  required  at 
an  interface  is  a  delicate  problem  that  remains  an  unsettled 
topic  in  the  literature.  From  the  physical  point  of  view, 
the  presence  of  a  pressure  in  the  electron  fluid  enables 
the  existence  of  longitudinal  waves ,  in  which  the  electric 
field  is  aligned  with  the  propagation  vector,  along  with  the 
conventional  transverse  waves.  That  is,  an  incident  field 
can  excite  two  independent  waves  inside  the  metal.  The 
well-known  Maxwell’s  boundary  conditions  are  not  sufficient 
to  uniquely  define  the  amplitudes  of  these  possible  waves.  In 
terms  of  a  hydrodynamic  description,  the  number  of  additional 
boundary  conditions  depends  on  the  given  equilibrium  charge 
density  profile  at  the  metal  boundaries.  For  our  problem,  we 
assume  the  equilibrium  electron  density  to  have  a  step-function 
profile  that  vanishes  outside  the  metal.  In  this  case,  only 
one  additional  boundary  condition  is  required  to  supplement 
Eqs.  (6)  to  obtain  a  unique  solution.82  Intuitively,  a  simple 
boundary  condition  may  be  ascertained  by  considering  the 
meaning  of  our  equations.  From  a  macroscopic  point  of  view, 
Maxwella  equations  require  the  normal  component  of  the 
vector  D  =  £q E  +  P  to  be  continuous  across  the  interface,  if 
no  external  charges  are  present,  while  the  normal  component 
of  the  electric  field  n  •  E  is  discontinuous  because  of  the 
induced  surface  charges  (n  is  a  unit  vector  normal  to  the 
boundary).  However,  the  nonlocal,  linear  terms  in  Eqs.  (6b) 
account  for  a  microscopic  description  of  the  charge  behavior 
at  the  surfaces  of  the  metal  region.  It  is  then  legitimate  to 
impose  the  continuity  of  the  normal  component  of  the  electric 
field  across  the  interface,  which  means  it  must  be  n  •  P  =  0 
at  the  metal-air  interface.  On  the  other  hand,  the  tangential 
component  of  the  polarization  vector,  iixP,  is,  in  general, 
nonzero.  More  formally,  the  same  result  may  be  derived  from 
the  continuity  equation  and  Gauss’s  theorem.  In  this  case,  one 
obtains  that  n  •  J  =  0  at  the  boundary.78,82 

Within  this  context  the  macroscopic  discontinuity  of 
n  •  E  translates  into  a  rapid  variation  of  the  microscopic 
fields  within  the  neighborhood  of  the  boundary  and  the 
term  V  •  P  will  give  indeed  the  distribution  of  the  induced 
charge  along  the  metal-air  transition.  Equations  (6)  coupled 
to  Eqs.  (8)  are  numerically  solved  using  the  finite-element 
method  implemented  in  the  commercial  software  COMSOL 
Multiphysics.84  The  weak  form  of  the  nonlocal  contribution 
derived  in  the  appendix  has  been  used. 

A.  Second-harmonic  generation  from  a  metal  nanowire 

To  illustrate  the  nature  of  the  various  nonlinear  processes 
that  derive  from  the  hydrodynamic  response  model,  we 
consider  a  p -polarized  plane  wave  incident  on  a  metal  wire  of 
circular  cross  section  as  depicted  in  Fig.  1(a).  It  is  interesting  to 
assess  the  influence  of  the  linear  pressure  term  on  the  induced 
charge  density.  One  would  expect  the  electron  distribution 
to  be  zero  in  the  bulk  region  and  rapidly  varying  near  the 
metal  surface,  where  it  should  reach  its  maximum  value. 
The  induced  charge  density  n  =  •  P  is  shown  in  Fig.  2, 

for  different  values  of  ft,  along  the  direction  normal  to  the 
metal-air  interface.  We  note  that  the  region  where  the  variation 
occurs  is  of  order  of  ~  0. 1  nm. 
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FIG.  1.  (Color  online)  (a)  A  plane  wave  of  wavelength  X  =  1064  nm  impinges  on  a  gold  wire  of  diameter  d  =  100  nm.  (b)  Absolute  value 
of  the  fundamental  electric  field.  The  second-harmonic  field  components  obtained  by  solving  Eqs.  (8)  and  (6)  are  shown  in  (c)  and  (d).  The 
following  values  were  used:  ra*  =  rae,  n0  =  5.7  x  1022  cm-3,  y  =  1.07  x  1014  s_1,  and  EF  =  5.5  eV. 


As  the  fundamental  wave  interacts  with  the  metal  wire, 
a  second-harmonic  field  given  by  the  nonlinear  polarization 
of  Eq.  (6b)  and  Eq.  (7)  is  generated  (see  Fig.  1).  It  is 
useful  to  closely  examine  the  contribution  of  the  various 
nonlinear  source  terms  to  the  generated  second  harmonic.  In 
Fig.  3  we  show  the  generated  second-harmonic  field  patterns 
corresponding  to  each  of  the  nonlinear  terms  in  Eq.  (6b), 
with  all  other  nonlinear  sources  turned  off.  The  corresponding 
second-harmonic  scattered  powers,  Qsh,  are  also  indicated  in 
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29.0  29.5  30.0  30.5  31.0 
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FIG.  2.  (Color  online)  Electron  density  profile  at  the  metal-air 
interface  for  different  values  of  /3  for  the  fundamental  field.  The 
geometry  was  discretized  by  a  triangular  mesh  with  a  0.025-nm 
maximum  element  size  at  the  metal  surface. 


Fig.  3.  If  all  the  nonlinear  sources  are  present  simultaneously 
the  total  converted  power  is  2sh  =  1-53  x  10“ 10  W,  a  value 
almost  twice  the  sum  given  by  each  of  the  individual 
contributions.  This  means  that  the  terms  interfere,  and  each 
term  may  establish  a  nonlinear  interaction  channel  that  could 
be  amplified  or  reduced  by  the  other  terms.  Consequently, 
it  is  not  easy  to  estimate  the  impact  of  an  individual  term 
on  the  generated  field,  which,  in  general,  will  depend  on  the 
geometrical  configuration  and  other  parameter  choices. 

It  is  also  useful  to  determine  the  relative  second  harmonic 
contributions  from  the  surface  and  the  bulk  of  the  geometry. 
With  respect  to  Eq.  (7),  surface  terms  may  be  identified  as  being 
proportional  to  V  -  Pi,  a  term  that  is  nonzero  only  near  the 
surface.  That  is,  the  Coulomb  (quadrupole-like)  and  pressure 
terms,  and  the  first  of  the  convective  terms,  may  be  taken  to 
be  purely  surface  terms,  while  the  Lorentz  force  is  a  purely 
bulk  term.  On  other  hand,  the  convective  term  proportional 
to  (Pi  •  V)Pi  contains  both  surface  and  bulk  contributions. 
In  Fig.  4  we  plot  the  absolute  values  of  the  transverse  and 
normal  component  of  the  vectors  (V  •  Pi)Pi  and  (Pi  •  V)Pi, 
respectively.  Near  the  surface  the  normal  components  of  both 
terms  behave  similarly,  while  only  the  latter  does  not  vanish  in 
the  bulk  region.  An  important  difference  can  be  observed  in  the 
transverse  component,  where  there  is  no  surface  contribution 
for  the  term  (Pi  •  V)Pi.  These  distinctions  qualitatively  hold 
for  different  types  of  nanoparticles.  However,  the  precise 
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FIG.  3.  (Color  online)  Second-harmonic  field  patterns  generated  by  the  different  nonlinear  sources  of  Eq.  (6b)  for  the  geometry  in  Fig.  1. 
The  total  scattered  second-harmonic  power  2Sh  is  shown  for  a  pump  peak  intensity  of  ~15  MW/cm2.  The  field  patterns  were  calculated  by 
considering  only  the  second-harmonic  source  term  of  interest  and  turning  off  all  others. 


fraction  associated  with  surface  or  bulk  contributions  depends 
on  the  specific  geometry  under  consideration. 

B.  Free-electron  limit 

The  hydrodynamic  description  of  the  electrons  inside  a 
metal  gives  a  fairly  accurate  description  of  linear  and  nonlinear 
processes  occurring  at  the  surface  of  metallic  structures. 
Though  microscopic  interactions  are  described  by  a  very 
simplistic  model,  the  simultaneous  manifestation  of  multiple 
scales  (the  angstrom  length  scale  of  the  electron-electron  inter¬ 
actions,  on  the  one  hand,  and  the  micron  scale  of  particle-field 
interaction,  on  the  other  hand)  makes  the  numerical  resolution 
of  the  electromagnetic  problem  quite  complex  and  ordinarily 
necessitates  considerable  computational  resources  even  for 
particles  whose  dimensions  are  a  few  tens  of  nanometers.  In 
Ref.  72  we  showed  that  it  is  possible  to  express  the  nonlinear 
surface  contribution  in  terms  of  the  bulk  polarization,  without 
having  to  solve  the  nonlocal  equations.  As  already  pointed  out 
by  Sipe  et  al .,35  if  the  electron  pressure  is  simply  dropped  off  by 
the  equations,  the  hydrodynamic  theory  of  second-harmonic 
generation  becomes  inherently  ambiguous.  This  is  because 
products  of  the  form  (V  •  P)E  are  not  well  defined  if  both  P 
and  E  are  discontinuous.72 

In  early  works  this  problem  was  circumvented  by  introduc¬ 
ing  phenomenological  coefficients  to  determine  the  weight  of 
the  different  nonlinear  contributions,34  or  through  an  effective 
plasma  frequency,35  that  incorporate  the  details  of  the  charge 
distribution  near  the  surface,  an  effect  that  was  neglected  in 
Ref.  70.  We  have  numerically  analyzed  the  effect  of  the  linear 
pressure  on  the  amount  of  converted  second-harmonic  energy. 


This  is  shown  in  Fig.  5  where  the  total  second-harmonic 
scattered  power  is  plotted  as  a  function  of  the  inverse  of  /3. 
Figure  5  shows  that  the  total  second-harmonic  generation 
converges  to  an  asymptotic  value  as  /3  tends  toward  zero. 
Typical  reported  values  for  are  on  the  order  of  106  m/s 
and  correspond  to  the  flat  region  of  the  curve,  where  the  actual 
converted  energy  does  not  differ  too  much  from  its  asymptotic 
value.  Exploring  the  limit  for  /3  0  seems,  then,  a  very 

good  way  to  get  an  approximate  solution  that  includes  the 
effects  of  nonlocality,  but  without  having  to  solve  the  complex 
nonlocal  equations.  Note  that  this  procedure  differs  totally 
from  assuming  /3  is  identically  zero. 


Normal  coordinate  (y)  [nm] 

FIG.  4.  (Color  online)  Comparison  between  the  convective  source 
terms  calculated  for  the  nanowire  of  Fig.  1.  The  nanowire  is  centered 
at  the  origin  and  the  fields  are  taken  along  the  direction  x  =  0  in  the 
neighborhood  of  the  metal-air  interface  (y  =  50  nm).  Note  that  the 
direction  of  interest  is  normal  to  the  metal  boundary. 
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FIG.  5.  (Color  online)  Total  second-harmonic  scattered  power  as 
a  function  of  the  inverse  of  the  parameter  /3  for  a  gold  nanowire  of 
diameter  d  =  30  nm. 


In  this  limit  the  surface  contributions  to  the  second- 
harmonic  polarization  may  be  approximated  by  an  effective 
nonlinear  current  sheet  at  the  surface  of  the  nanoparticle 
given  by72 
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where  the  unit  vectors  n  and  t  are  normal  and  parallel  to  the 
metal  interface,  respectively,  and  =  n  Pi  and  =  t  •  Pi. 
These  currents  are  related  to  the  polarization  values  in  the  bulk 
region  and  do  not  require  the  resolution  of  nonlocal  equations. 
This  approach  provides  a  good  description  of  the  second- 
harmonic  generation  process  that  may  be  easily  implemented 
for  full  three-dimensional  (3D)  simulations.  However,  since 
in  this  local  description  of  surface  nonlinear  contributions,  the 
derivatives  in  the  direction  parallel  to  the  surface  are  neglected, 
Eq.  (9)  is  not  expected  to  hold  if,  for  example,  one  considers 
geometries  with  sharp  asperities,  in  which  field  derivatives 
can  be  very  important.72  More  generally,  the  approximation 
is  valid  as  long  as  nonlocal  effects  do  not  change  the  linear 
response  of  the  system. 

As  a  test  case,  we  consider  second-harmonic  generation 
from  a  thick  metal  film.  In  Ref.  84  experimental  data  for  a  pump 
field  tuned  at  1064  nm  and  incident  on  a  400-nm  silver  film 
are  reported  for  three  different  incident/detected  polarization 
configurations:  (i)  TM-incident  and  TM-detected  polarization; 
(ii)  TE-incident  and  TM-detected  polarization,  and  (iii)  45° - 
incident  and  TE-detected  polarization.  In  Fig.  6  we  compared 
the  results  obtained  using  Eq.  (9)  with  those  obtained  through 
the  full  resolution  of  Eqs.  (6).  For  completeness  the  experi¬ 
mental  data  are  also  reported.  Both  simulation  results  agree 
qualitatively  and  quantitatively  very  well  in  all  cases.  A  slight 
difference  between  the  two  numerical  results  may  be  observed 
for  the  case  of  TM-incident  and  TM-detected  polarization 
[see  Fig.  6(a)].  This  slight  difference  arises  because  in  that 
configuration  the  impact  of  the  electron  pressure  is  more 
significant.  However,  the  qualitative  agreement  is  still  very 
good. 


III.  THREE-DIMENSIONAL  NANOPARTICLES 

In  the  previous  paragraph  we  showed  how  second- 
harmonic  generation  may  be  handled  without  solving  the 
entire  electromagnetic-hydrodynamic  problem.  That  is,  the 


10-i4  TE/TM  incident,  TE  detected 


FIG.  6.  (Color  online)  Second-harmonic  conversion  efficiency 
for  a  silver  film  as  a  function  of  the  incident  angle  6  as  depicted  in 
the  inset,  under  various  pumping/detection  polarization  conditions. 
The  fundamental  field  is  tuned  at  the  wavelength  X  =  1064  nm. 
Experimental  data  are  taken  from  Ref.  84,  where  the  average  peak 
intensity  is  estimated  about  80  MW/cm2.  The  shapes  of  both 
theoretical  curves  agree  well  with  the  experimental  data  in  all  cases,  if 
we  choose  the  following  incident  peak  intensities:  (a)  ~  85  MW /cm2, 
(b)  ~  60  MW/cm2,  and  (c)  ~  110  MW/cm2.  The  effective  electron 
mass  has  been  fixed  to  m*  =  0.65me  and  n0  =  3.5  x  1022  cm-3, 
y  =  4.6  x  1013  s-\  and  E¥  =  5.5  eV. 

volume  sources  are  calculated  by  assuming  identically  zero, 
while  the  surface  nonlinear  currents  are  given  by  Eq.  (9). 
In  this  section  we  will  use  this  approach  to  study  3D  gold 
nanoparticles.  For  the  sake  of  simplicity,  we  assume  the 
nanoparticles  are  surrounded  by  air.  The  structure  studied 
extends  periodically  in  the  x  and  y  directions,  as  depicted 
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FIG.  7.  (Color  online)  Second-harmonic  conversion  efficiency  for  different  3D  gold  nanoparticles.  All  the  particles  are  20  nm  thick  and 
the  geometrical  parameters  have  been  chosen  to  have  a  resonance  around  A.FF  =  1.5  /zm,  where  the  fundamental  field  is  tuned.  The  linear 
transmissions  at  normal  incidence  (solid  line)  for  the  variety  of  nanoparticles  are  shown  in  the  second  column.  In  (a)  and  (b)  the  transmission  at 
oblique  incidence  (dashed  line)  has  been  calculated  for  6  =  40  and  0  =  10,  respectively.  The  transmission  for  an  incident  electric  field  polarized 
along  the  y  direction  is  also  shown  for  the  U-shaped  nanoparticle  (dot-dash  line).  The  vertical  dotted  lines  indicate  the  fundamental  and  the 
second-harmonic  wavelength,  respectively.  The  second-harmonic  conversion  efficiency  as  a  function  of  the  incident  angle  0  has  been  calculated 
using  Eq.  (9)  to  take  into  account  the  nonlinear  surface  contribution.  The  pumping  electric  field  is  p  polarized  (Ey  =  0)  with  a  peak  intensity 
of  ~55  MW/cm2.  The  following  values  for  the  parameters  have  been  used:  m*  =  me,  n0  =  5.7  x  1022  cm-3,  and  y  =  1.07  x  1014  s-1. 


in  the  inset  of  Fig.  7,  so  only  a  single  unit  cell  is  needed  in 
the  computational  space.  To  avoid  possible  numerical  artifacts 
due  to  the  field  localization  near  metal  corners,  we  considered 
rounded  comer  geometries  with  a  radius  of  curvature  of  5  nm. 
The  geometrical  parameters  were  chosen  so  the  nanoparticles 
would  display  a  resonance  around  A.FF  =  1.5  /z m,  where  the 
fundamental  field  is  tuned. 


Our  results  are  summarized  in  Fig.  7  for  the  different  kinds 
of  nanoparticles.  The  conversion  efficiencies  assume  an  aver¬ 
age  pump  intensity  of  ~55  MW / cm2  with  the  electric  pumping 
field  lying  on  the  xz  plane.  We  find  qualitatively  good  agree¬ 
ment  with  the  experimental  data  of  Ref.  58,  where  second-  and 
third-harmonic  generation  were  experimentally  investigated 
for  a  variety  of  gold  nanoparticles.  The  authors  considered 
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arrays  of  unit  cells  ranging  from  300  to  600  nm  constituted  of 
gold  nanoparticles  deposited  on  a  thick  glass  substrate  under 
normal  pump  incidence.  Even  if  a  quantitative  comparison 
between  our  simulations  and  the  experimental  data  cannot 
be  made,  the  relative  efficiencies  normalized  with  respect  to 
the  U-shaped  nanoparticle  efficiency  are,  nevertheless,  in  very 
good  agreement.  For  the  U-shaped  nanoparticle  at  normal 
incidence,  we  find  a  conversion  efficiency  of  ~6.9  x  10-10. 
This  should  be  compared  to  ~2.0  x  10“ 11  of  the  experimental 
data.58,85  The  T-shaped  nanoparticle  shows  an  efficiency  of 
~1.1  x  10-11  for  0  =  0°,  or  about  the  1.7%  with  respect  of 
the  U-shaped  nanoparticle  (2.1%  for  the  experimental  data), 
and  the  I -shaped  nanoparticle  has  a  ~0  conversion  efficiency, 
as  one  might  expect  due  to  the  absence  of  symmetry  breaking. 
It  is  interesting  to  observe  the  second-harmonic  conversion 
efficiency  as  a  function  of  pumping  angle  6  for  the  three 
kinds  of  nanoparticles.  The  geometry  is  depicted  in  the  third 
column  of  Fig.  7.  While  the  U-shaped  nanoparticle  has  its 
conversion  peak  at  0°,  both  the  I-  and  T-shaped  nanoparticles 
show  maximum  efficiency  at  oblique  incidence.  In  particular, 
the  T -shaped  nanoparticle  qualitatively  behaves  as  a  bare  metal 
film,  reaching  a  maximum  second-harmonic  enhancement  at 
large  angles  (i 9  ~  70°),  as  shown  in  Fig.  7(b).  The  l-shaped 
nanoparticle  shows  a  more  interesting  response.  It  reaches  a 
conversion  efficiency  that  overtakes  the  enhancement  obtained 
with  the  U-shaped  nanoparticle  around  6  —  40° .  An  analysis  of 
the  linear  transmission  shows  that  the  second  resonance  [see 


FIG.  8.  (Color  online)  (a)  Fundamental  and  second-harmonic 
field  distributions  for  150-nm-wide  and  20-nm- thick  gold  nanocres¬ 
cent  arranged  in  a  periodical  array,  (b)  Electric  field  polarization  state 
for  the  fundamental  and  second-harmonic  fields,  respectively.  The 
polarization  state  of  the  generated  field  is  flipped  with  respect  to  the 
pumping  field. 


Fig.  7(a)]  redshifts  as  6  increases,  attaining  the  wavelength 
A.ff/2  for  0  ~  40°,  exactly  where  the  maximum  conversion 
occurs.  The  energy  flowing  from  the  pump  to  the  second- 
harmonic  field  is  then  catalyzed,  thanks  to  the  presence  of  this 
double  resonance,  as  shown  in  Fig.  7(a).  The  situation  differs 
slightly  for  the  T-shaped  nanoparticle,  where  the  resonances 
are  matched  only  for  wide  angles,  when  the  v  component  of 
the  electric  field,  the  one  that  drives  the  resonance,  is  sensibly 
reduced. 

Another  interesting  structure  to  study  is  the  nanocrescent. 
This  structure  has  been  shown  to  be  able  to  efficiently  harvest 
light  over  a  broadband  spectrum.86  Moreover,  its  inherent 
asymmetry  makes  it  a  good  candidate  for  second-harmonic 
generation  enhancement.  In  Fig.  8(a)  we  show  the  funda¬ 
mental  and  second-harmonic  field  distributions  for  120-nm- 
wide  periodically  arranged  nanocrescents.  The  nanocrescent 
was  designed  to  reach  a  maximum  sharpness  defined  by  a 
minimum  allowed  size  8  and  to  have  a  gap  of  size  28.  The 
structure  is  excited  by  a  fundamental  wave  polarized  along 
the  horizontal  direction,  as  shown  in  Fig.  8(b).  Similarly  to 
the  U-shaped  nanoparticle,  the  symmetry  breaking  occurs 
along  the  y  direction,  the  polarization  state  of  the  generated 
field  results  then  rotated  by  90°.  The  structure  shows  a 
maximum  conversion  efficiency  at  normal  incidence  that  can 
enhance  the  second-harmonic  generation  process  up  to  two 
orders  of  magnitude  over  that  of  the  U-shaped  nanoparticle 
efficiency.  In  particular,  we  obtained  conversion  efficiencies 
of  r]  ~  6.0  x  10-8,  rj  ~  3.5  x  10-8,  and  r\  ~  5.1  x  10-9  for 
8  =  2,3,5  nm,  respectively,  corresponding  to  a  pumping  peak 
intensity  of  ~55  MW/cm2. 


IV.  CONCLUSION 

We  have  studied  the  second-harmonic  generation  from 
metal-based  nanostructures.  The  nonlinear  optical  response 
of  metal  was  described  using  the  hydrodynamic  model, 
which  introduced  a  source  of  nonlocality  in  the  dispersion 
relation,  due  to  electron  gas  pressure.  We  have  discussed 
its  influence  on  second-harmonic  generation  and  numerically 
showed  that,  as  the  pressure  term  tends  to  zero,  the  amount  of 
converted  second-harmonic  field  tends  to  an  asymptotic  value, 
leading  to  the  possibility  of  expressing  the  nonlinear  surface 
contributions  as  a  function  of  the  values  of  the  polarization 
vector  in  the  bulk  region.  This  development  simplifies  the 
investigation  of  second-harmonic  generation  process  in  full 
3D  metal  structures. 

Recent  experimental  work  has  shown  significant  enhance¬ 
ment  of  the  electric  field  in  plasmonic  systems  when  metal 
nanoparticles  are  strongly  coupled  to  a  metal  film.76,87  As  the 
distances  between  the  nanoparicles  and  the  film  approach  the 
scale  of  fraction  of  a  nanometer  nonlocal  effects  are  no  longer 
negligible,  and  nonlinear  effects  become  stronger.  We  believe 
that  the  ability  to  solve  the  full  electromagnetic-hydrodynamic 
problem  is  crucial  for  this  kind  of  system  and  will  be  the  subject 
of  our  future  investigations. 
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APPENDIX:  NONLOCAL  WEAK  CONTRIBUTION 

Consider  a  differential  problem  of  the  form  Lu  =  /  with  L 
a  linear  differential  operator  and  /  an  arbitrary  function.  The 
associated  weak  form  is  then 


/  cpLudV  =  f  (pfdV, 

Jn  Jn 


(Al) 


with  u\dQ  =  uo,  where  tp  is  en  element  of  a  set  of  arbitrary 
functions  called  test  functions.  Generally,  the  form  of  the 
integrals  of  Eq.  (Al)  allows  us  to  decrease  the  order  of 
derivatives  of  the  operator  L.  Here,  we  present  a  suitable  weak 
form  for  the  problem  of  Eqs.  (6). 

For  the  sake  of  simplicity,  we  rewrite  here  only  differential 
higher-order  terms,  since  lower  ones  remain  unchanged. 
Consider,  first,  the  term  associated  with  the  linear  electron 
pressure  V  (V  •  P).  Multiplying  for  the  vectorial  test  function 
O  and  integrating  over  the  volume  as  in  Eq.  (Al),  we  obtain 


/ 

Jn 


O  .  [V(V  .  P )]dV. 


(A2) 


Integrating  by  parts  and  using  the  divergence  theorem,  we 
obtain 


/ 

Jn 


O  •  [V(V  .  P )]dV 


=  [  (V  .  P)n  •  O  dA  -  f  (V  .  cD)(V  .  P )dv.  (A3) 

J  3  £2  J 

Finally,  considering  that  O  •  n  =  0  on  the  boundaries,  we 
obtain 

f  o  .  [V(V  .  P )]dV  =  -  [  (V  •  0)(V  .  P )dV.  (A4) 

J  n  J  n 

Analogously,  we  obtain,  for  the  nonlinear  pressure  term, 
(V  •  P)  V  (V  •  P),  of  Eq.  (6b)  the  following  result: 


f  <D  •  [(V  •  P)V(V  •  P )]dV 
Jn 


where  we  used  the  identity 


-\L 


(V  •  4>)(V  -P fdV, 


(V  -P)V(V  -P)  =  ±V(V  -P)2 


(A5) 


(A6) 


The  surface  currents  given  by  Eq.  (9)  have  been  imple¬ 
mented  by  using  surface  contributions. 


*  cri  stian .  cir  aci  @  duke .  edu 

lD.  Schurig,  J.  J.  Mock,  B.  J.  Justice,  S.  A.  Cummer,  J.  B.  Pendry, 
A.  F.  Starr,  and  D.  R.  Smith,  Science  314,  977  (2005). 

2D.  R.  Smith  and  N.  Kroll,  Phys.  Rev.  Lett.  85,  2933  (2000). 

3I.  V.  Shadrivov,  A.  B.  Kozyrev,  D.  van  der  Weide,  and  Y.  S.  Kivshar, 
Opt.  Express  16,  20266  (2008). 

4A.  A.  Zharov,  I.  V.  Shadrivov,  and  Y.  S.  Kivshar,  Phys.  Rev.  Lett. 
91,  037401  (2003). 

5S.  O’Brien,  D.  Mcpeake,  S.  A.  Ramakrishna,  and  J.  B.  Pendry, 
Phys.  Rev.  B  69,  241101  (2004). 

6E.  Poutrina,  D.  Huang,  and  D.  R.  Smith,  New  J.  Phys.  12,  093010 
(2010). 

7S.  Larouche  and  D.  R.  Smith,  Opt.  Commun.  283,  1621  (2010). 

8 A.  Rose,  S.  Larouche,  D.  Huang,  E.  Poutrina,  and  D.  R.  Smith, 
Phys.  Rev.  E  82,  036608  (2010). 

9J.  Zhou,  T.  Koschny,  M.  Kafesaki,  E.  N.  Economou,  J.  B.  Pendry, 
and  C.  M.  Soukoulis,  Phys.  Rev.  Lett.  95,  223902  (2005). 

10R.  W.  Boyd  and  G.  L.  Fischer,  in  Encyclopedia  of  Materials: 
Science  and  Technology  ( Elsevier  Science  Ltd.,  Amsterdam,  2001), 
pp.  6237-6244. 

nP.  A.  Franken,  G.  Weinreich,  C.  W.  Peters,  and  A.  E.  Hill,  Phys. 
Rev.  Lett.  7,  118  (1960). 

12N.  Bloembergen  and  P.  S.  Pershan,  Phys.  Rev.  128,  606  (1961). 

13P.  S.  Pershan,  Phys.  Rev.  130,  919  (1963). 

14P.  A.  Franken  and  J.  Ward,  Rev.  Mod.  Phys.  35,  23  (1962). 

15N.  Bloembergen,  Proc.  IEEE  51,  124  (1963). 

16E.  Adler,  Phys.  Rev.  134,  A728  (1964). 

17N.  Bloembergen,  R.  Chang,  S.  Jha,  and  C.  Lee,  Phys.  Rev.  174,  813 
(1968). 

18F.  Brown,  R.  Parks,  and  A.  Sleeper,  Phys.  Rev.  Lett.  14, 1029  (1965). 


19S.  Jha,  Phys.  Rev.  Lett.  15,  412  (1965). 

20H.  Sonnenberg  and  H.  Heffner,  J.  Opt.  Soc.  Am.  A  58,  209 
(1968). 

21  J.  C.  Quail  and  H.  J.  Simon,  Phys.  Rev.  B  31,  4900  (1985). 

22 J.  Coutaz,  D.  Maystre,  M.  Neviere,  and  R.  Reinisch,  J.  Appl.  Phys. 
62,  1529  (1987). 

23H.  J.  Simon,  C.  Huang,  J.  C.  Quail,  and  Z.  Chen,  Phys.  Rev.  B  38, 
7408  (1988). 

24D.  Krause,  C.  W.  Teplin,  and  C.  T.  Rogers,  J.  Appl.  Phys.  96,  3626 
(2004). 

25S.  Jha,  Phys.  Rev.  140,  A2020  (1965). 

26N.  Bloembergen  and  Y.  Shen,  Phys.  Rev.  141,  298  (1966). 

27P.  Guyot-Sionnest,  W.  Chen,  and  Y.  R.  Shen,  Phys.  Rev.  B  33,  8254 
(1986). 

28F.  Brown  and  R.  Parks,  Phys.  Rev.  Lett.  16,  507  (1966). 

29 A.  Liebsch,  Phys.  Rev.  Lett.  61,  1233  (1988). 

30 J.  I.  Dadap,  J.  Shan,  K.  B.  Eisenthal,  and  T.  F.  Heinz,  Phys.  Rev. 
Lett.  83,  4045  (1999). 

31M.  I.  Stockman,  D.  J.  Bergman,  C.  Anceau,  S.  Brasselet,  and 
J.  Zyss,  Phys.  Rev.  Lett.  92,  057402  (2004). 

32K.  Li,  M.  I.  Stockman,  and  D.  J.  Bergman,  Phys.  Rev.  B  72,  153401 
(2005). 

33G.  Bachelier,  I.  Russier- Antoine,  E.  Benichou,  C.  Jonin,  and  P.-F. 

Brevet,  J.  Opt.  Soc.  Am.  B  25,  955  (2008). 

34J.  Rudnick  and  E.  Stern,  Phys.  Rev.  B  4,  4274  (1971). 

35 J.  E.  Sipe,  V.  C.  Y.  So,  M.  Fukui,  and  G.  I.  Stegeman,  Phys.  Rev.  B 
21,  4389  (1980). 

36M.  Corvi  and  W.  L.  Schaich,  Phys.  Rev.  B  33,  3688  (1986). 

37C.  D.  Hu,  Phys.  Rev.  B  40,  7520  (1989). 

38M.  Airola,  Y.  Liu,  and  S.  Blair,  J.  Opt.  A  7,  SI  18  (2005). 


115451-9 


CIRACI,  POUTRINA,  SCALORA,  AND  SMITH 


PHYSICAL  REVIEW  B  86,  115451  (2012) 


39 A.  Lesuffleur,  L.  K.  S.  Kumar,  and  R.  Gordon,  Appl.  Phys.  Lett. 
88,  261104  (2005). 

40  J.  A.  H.  van  Nieuwstadt,  M.  Sandtke,  R.  H.  Harmsen,  F.  B.  Segerink, 
J.  C.  Prangsma,  S.  Enoch,  and  L.  Kuipers,  Phys.  Rev.  Lett.  97, 
146102  (2006). 

41M.  A.  Vincenti,  V.  Petruzzelli,  A.  D’Orazio,  F.  Prudenzano,  M.  J. 
Bloemer,  N.  Akozbek,  and  M.  Scalora,  J.  Nanophotonics  2, 021851 
(2008). 

42T.  Xu,  X.  Jiao,  and  S.  Blair,  Opt.  Express  17,  23582  (2009). 

43M.  A.  Vincenti,  D.  de  Ceglia,  V.  Roppo,  and  M.  Scalora,  Opt. 
Express  19,  2067  (2011). 

44M.  C.  Larciprete,  A.  Belardini,  M.  G.  Cappeddu,  D.  de  Ceglia, 
M.  Centini,  E.  Fazio,  C.  Sibilia,  M.  J.  Bloemer,  and  M.  Scalora, 
Phys.  Rev.  A  77,  013809  (2008). 

45G.  D’Aguanno,  M.  C.  Larciprete,  N.  Mattiucci,  A.  Belardini,  M.  J. 
Bloemer,  E.  Fazio,  O.  Buganov,  M.  Centini,  and  C.  Sibilia,  Phys. 
Rev.  A  81,  013834  (2010). 

46K.  D.  Ko,  A.  Kumar,  K.  H.  Fung,  R.  Ambekar,  G.  L.  Liu,  N.  Fang, 
and  K.  C.  Toussaint,  Nano  Lett.  11,  61  (2011). 

47 A.  Benedetti,  M.  Centini,  C.  Sibilia,  and  M.  Bertolotti,  J.  Opt.  Soc. 
Am.  B  27,  408  (2010). 

48 A.  Belardini,  M.  C.  Larciprete,  M.  Centini,  E.  Fazio,  C.  Sibilia, 
M.  Bertolotti,  A.  Toma,  D.  Chiappe,  and  FB.de  Mongeot,  Opt. 
Express  17,  3603  (2009). 

49J.  Renger,  R.  Quidant,  N.  van  Hulst,  and  L.  Novotny,  Phys.  Rev. 
Lett.  104,  046803  (2010). 

50M.  W.  Klein,  C.  Enkrich,  M.  Wegener,  and  S.  Linden,  Science  313, 
502  (2006). 

51F.  B.  P.  Niesler,  N.  Feth,  S.  Linden,  J.  Niegemann,  J. 
Gieseler,  K.  Busch,  and  M.  Wegener,  Opt.  Lett.  34,  1997 
(2009). 

52C.  C.  Neacsu,  G.  A.  Reider,  and  M.  B.  Raschke,  Phys.  Rev.  B  71, 
201402(R)  (2005). 

53J.  Nappa,  G.  Revillod,  I.  Russier-Antoine,  E.  Benichou,  C.  Jonin, 
and  P.  F.  Brevet,  Phys.  Rev.  B  71,  165407  (2005). 

54J.  Shan,  J.  I.  Dadap,  I.  Stiopkin,  G.  A.  Reider,  and  T.  F.  Heinz,  Phys. 
Rev.  A  73,  023819  (2006). 

55M.  Zavelani-Rossi,  M.  Celebrano,  P.  Biagioni,  D.  Polli,  M.  Finazzi, 
L.  Duo,  G.  Cerullo,  M.  Labardi,  M.  Allegrini,  J.  Grand,  and 
P.  Adam,  Appl.  Phys.  Lett.  92,  093119  (2008). 

56  J.  Butet,  J.  Duboisset,  G.  Bachelier,  I.  Russier-Antoine, 
E.  Benichou,  C.  Jonin,  andP.-F.  Brevet,  Nano  Lett.  10, 1717  (2010). 
57B.  K.  Canfield,  S.  Kujala,  K.  Jefimovs,  J.  Turunen,  andM.  Kauranen, 
Opt.  Express  12,  5418  (2004). 

58M.  W.  Klein,  M.  Wegener,  N.  Feth,  and  S.  Linden,  Opt.  Express  15, 
5238  (2006). 

59M.  D.  McMahon,  R.  Lopez,  R.  F.  Haglund,  E.  A.  Ray,  and  P.  H. 

Bunton,  Phys.  Rev.  B  73,  041401(R)  (2006). 

60B.  K.  Canfield,  H.  Husu,  J.  Laukkanen,  B.  Bai,  M.  Kuittinen, 
J.  Turunen,  and  M.  Kauranen,  Nano  Lett.  7,  1251  (2007). 

61 C.  Hubert,  L.  Billot,  P.  Adam,  R.  Bachelot,  P.  Royer,  J.  Grand, 
D.  Gindre,  K.  D.  Dorkenoo,  and  A.  Fort,  Appl.  Phys.  Lett.  90, 
181105  (2007). 


62N.  Feth,  S.  Linden,  M.  W.  Klein,  M.  Decker,  F.  B.  P.  Niesler, 
Y.  Zeng,  W.  Hoyer,  J.  Liu,  S.  W.  Koch,  J.  V.  Moloney,  and 
M.  Wagener,  Opt.  Lett.  33,  1975  (2008). 

63H.  Husu,  B.  K.  Canfield,  J.  Laukkanen,  B.  Bai,  M.  Kuittinen, 
J.  Turunen,  and  M.  Kauranen,  Appl.  Phys.  Lett.  93, 183115  (2008). 

64E.  Kim,  F.  Wang,  W.  Wu,  Z.  Yu,  and  Y.  R.  Shen,  Phys.  Rev.  B  78, 
113102  (2008). 

65 S.  Kujala,  B.  K.  Canfield,  M.  Kauranen,  Y.  Svirko,  and  J.  Turunen, 
Opt.  Express  16,  17196  (2008). 

66 V.  K.  Valev,  N.  Smisdom,  A.  V.  Silhanek,  B.  de  Clercq,  W.  Gillijns, 

M.  Ameloot,  V.  V.  Moshchalkov,  and  T.  Verbiest,  Nano  Lett.  9, 
3945  (2009). 

67  V.  K.  Valev,  A.  V.  Silhanek,  N.  Smisdom,  B.  de  Clercq,  W.  Gillijns, 
O.  Aktsipetrov,  M.  Ameloot,  V.  V.  Moshchalkov,  and  T.  Verbiest, 
Opt.  Express  18,  8286  (2010). 

68V.  K.  Valev,  A.  V.  Silhanek,  N.  Verellen,  W.  Gillijns,  P.  van  Dorpe, 
O.  A.  Aktsipetrov,  G.  A.  E.  Vandenbosch,  V.  V.  Moshchalkov,  and 
T.  Verbiest,  Phys.  Rev.  Lett.  104,  127401  (2010). 

69T.  Utikal,  T.  Zentgraf,  T.  Paul,  C.  Rockstuhl,  F.  Lederer,  M.  Lippitz, 
and  H.  Giessen,  Phys.  Rev.  Lett.  106,  133901  (2011). 

70Y.  Zeng,  W.  Hoyer,  J.  Liu,  S.  W.  Koch,  and  J.  V.  Moloney,  Phys. 
Rev.  B  79,  235109  (2009). 

71 M.  Scalora,  M.  A.  Vincenti,  D.  de  Ceglia,  V.  Roppo,  M.  Centini, 

N.  Akozbek,  and  M.  J.  Bloemer,  Phys.  Rev.  A  82,  043828 

(2010). 

72C.  Ciraci,  E.  Poutrina,  M.  Scalora,  and  D.  R.  Smith,  Phys.  Rev.  B 
85,  201 403 (R)  (2012). 

73 A.  Eguiluz  and  J.  Quinn,  Phys.  Rev.  B  14,  1347  (1976). 

74N.  Crouseilles,  P.-A.  Hervieux,  and  G.  Manfredi,  Phys.  Rev.  B  78, 
155412  (2008). 

75R.  W.  Boyd,  Nonlinear  Optics  (Academic  Press,  San  Diego,  CA, 
2006). 

76C.  Ciraci,  R.  T.  Hill,  J.  J.  Mock,  Y.  A.  Urzhumov,  A.  I.  Fernandez - 
Dommguez,  S.  A.  Maier,  J.  B.  Pendry,  A.  Chilkoti,  and  D.  R.  Smith, 
Science  337,  1072  (2012). 

77R.  Ruppin,  Opt.  Commun.  190,  205  (2001). 

78 S.  Raza,  G.  Toscano,  A.-P.  Jauho,  M.  Wubs,  and  N.  A.  Mortensen, 
Phys.  Rev.  B  84,  121412  (2011). 

79 J.  M.  Mcmahon,  S.  K.  Gray,  and  G.  C.  Schatz,  Phys.  Rev.  Lett.  103, 
097403  (2009). 

80P.  Halevi  and  R.  Fuchs,  J.  Phys.  C  17,  3889  (1984). 

81D.  Maystre,  M.  Neviere,  and  R.  Reinisch,  Appl.  Phys.  A  39,  115 
(1986). 

82P.  Jewsbury,  J.  Phys.  F  11,  195  (1981). 

83  COMSOL  Multiphysics,  http://www.comsol.com/. 

84K.  O’Donnell  and  R.  Torre,  New  J.  Phys.  7,  154  (2005). 

85M.  W.  Klein,  M.  Wegener,  N.  Feth,  and  S.  Linden,  Opt.  Express  16, 
8055  (2008). 

86 A.  Aubry,  D.  Y.  Lei,  A.  I.  Fernandez-Dommguez,  Y.  Sonnefraud, 
S.  A.  Maier,  and  J.  B.  Pendry,  Nano  Lett.  10,  2574  (2010). 

87R.  T.  Hill,  J.  J.  Mock,  Y.  A.  Urzhumov,  D.  S.  Sebba,  S.  J.  Oldenburg, 
S.-Y.  Chen,  A.  A.  Lazarides,  A.  Chilkoti,  and  D.  R.  Smith,  Nano 
Lett.  10,  4150  (2010). 


115451-10 


